function out = getEquityPrem(model,modelRealYields,outSim,T,consClaim)

%outSim  = simProjection(model,setupFilter.shocks);
posMuz  = find(strcmp(model.labelx,'$\mu _{z,t}$')); 
if consClaim == 1
    posDiv  = find(strcmp(model.labely,'$c_t$')); 
    posEq   = find(strcmp(model.labely,'$pc^m_t$')); 
else
    posDiv  = find(strcmp(model.labely,'$div_t$')); 
    posEq   = find(strcmp(model.labely,'$p^m_t$')); 
end

logmuz_cu = outSim.x_cu(posMuz,:)+log(model.params.MUZss);
div_cu  = outSim.y_cu(posDiv,:);
eq_cu   = outSim.y_cu(posEq,:);
% Equity return 
retEq_cu= log([1 (exp(div_cu(2:end))+exp(eq_cu(2:end)))./exp(eq_cu(1:end-1))])+logmuz_cu;
% Real yield
if model.pruningOn == 1
    realYields_cu = getYields(modelRealYields,outSim.xAll_cu);
    yields_cu     = getYields(model,outSim.xAll_cu);
elseif model.pruningOn == 0
    realYields_cu = getYields(modelRealYields,outSim.x_cu);
    yields_cu     = getYields(model,outSim.x_cu);
end

% Equity premium
levFactor       = 1;
exReturn_cu     = levFactor*retEq_cu-[realYields_cu(1,1) realYields_cu(1,1:end-1)];
out.meanEqPrem  = mean(exReturn_cu*4);
out.stdEqPrem   = std(exReturn_cu)*sqrt(4);
out.meanRetEq   = mean(levFactor*retEq_cu*4);
out.stdRetEq    = std(levFactor*retEq_cu)*sqrt(4);

% Computing annual price-dividend ratio
numSim          = size(exReturn_cu,2);
logz_cu         = cumsum(logmuz_cu-log(model.params.MUZss));
% Cannot include the deterministic trend because we then get overflow. But
% the trend should cancel out for both prices and dividends.
% The leverage factor is 2 i.e. div^2
div_quarterly   = exp(div_cu).*exp(logz_cu);
div_annual      = zeros(1,numSim);
exReturn_annual = zeros(1,numSim);
for s=1:numSim
    if s == 1
        div_annual(1,s)      = div_quarterly(1,s)*4;
        exReturn_annual(1,s) = exReturn_cu(1,s)*4;
    elseif s == 2
        div_annual(1,s)      = (div_quarterly(1,s)+div_quarterly(1,s-1))*2;
        exReturn_annual(1,s)   = (exReturn_cu(1,s)+exReturn_cu(1,s-1))*2;
    elseif s == 3
        div_annual(1,s) = div_quarterly(1,s)+div_quarterly(1,s-1)...
            +div_quarterly(1,s-2)*4/3;
        exReturn_annual(1,s) = (exReturn_cu(1,s)+exReturn_cu(1,s-1)+exReturn_cu(1,s-2))*4/3;
    elseif s>3
        div_annual(1,s) = div_quarterly(1,s)+div_quarterly(1,s-1)...
            +div_quarterly(1,s-2)+div_quarterly(1,s-3);
        exReturn_annual(1,s) = (exReturn_cu(1,s)+exReturn_cu(1,s-1)+exReturn_cu(1,s-2)+exReturn_cu(1,s-3));
    end
end
pd              = exp(eq_cu).*exp(logz_cu)./div_annual;    %the deterministic trend cancel approximately out
out.meanLogPd   = mean(log(pd));
out.stdLogPd    = std(log(pd));
if numSim < 1000
    out.pd      = log(pd);
    out.rmReal  = retEq_cu;
    out.exReturn_cu = exReturn_annual;
end

% Predicting excess returns by the slope of the yield curve (and pd)
maxHorizon         = 12;
regExR_on_pd       = zeros(maxHorizon,2);
regExR_on_slope    = zeros(maxHorizon,2);
regExR_on_slope_pd = zeros(maxHorizon,3);
%exR_z              = (exReturn_cu'-mean(exReturn_cu))./std(exReturn_cu);
exR_z              = (exReturn_annual'-mean(exReturn_annual))./std(exReturn_annual);
log_pd_z           = (log(pd)'-mean(log(pd)))/std(log(pd));
slope              = (yields_cu(end,:)-yields_cu(4,:))*4;
slope_z            = (slope'-mean(slope))/std(slope);
for j=1:maxHorizon
    results = predictRegress(exR_z,log_pd_z,j,j*2);
    regExR_on_pd(j,:)   = [results.beta(2,1) results.rsqr];
    results = predictRegress(exR_z,slope_z,j,j*2);
    regExR_on_slope(j,:)   = [results.beta(2,1) results.rsqr];
    results = predictRegress(exR_z,[slope_z log_pd_z],j,j*2);
    regExR_on_slope_pd(j,:)   = [results.beta(2,1) results.beta(3,1) results.rsqr];
end
out.regExR_on_pd = regExR_on_pd;
out.regExR_on_slope = regExR_on_slope;
out.regExR_on_slope_pd = regExR_on_slope_pd;

%% accounting for finite sample biases
numRep = floor(numSim/T);
regExR_on_pd       = zeros(maxHorizon,2,numRep);
regExR_on_slope    = zeros(maxHorizon,2,numRep);
regExR_on_slope_pd = zeros(maxHorizon,3,numRep);
for s=1:numRep
    exR_z              = (exReturn_annual(1+(s-1)*T:s*T)'-mean(exReturn_annual(1+(s-1)*T:s*T)))./std(exReturn_annual(1+(s-1)*T:s*T));
    log_pd_z           = (log(pd(1+(s-1)*T:s*T))'-mean(log(pd(1+(s-1)*T:s*T))))/std(log(pd(1+(s-1)*T:s*T)));
    slope_z            = (slope(1+(s-1)*T:s*T)'-mean(slope(1+(s-1)*T:s*T)))/std(slope(1+(s-1)*T:s*T));
    for j=1:maxHorizon
        results = predictRegress(exR_z,log_pd_z,j,0);
        regExR_on_pd(j,:,s)   = [results.beta(2,1) results.rsqr];
        results = predictRegress(exR_z,slope_z,j,0);
        regExR_on_slope(j,:,s)   = [results.beta(2,1) results.rsqr];
        results = predictRegress(exR_z,[slope_z log_pd_z],j,0);
        regExR_on_slope_pd(j,:,s)   = [results.beta(2,1) results.beta(3,1) results.rsqr];
    end
end
out.regExR_on_pd_finiteSample = mean(regExR_on_pd,3);
out.regExR_on_slope_finiteSample = mean(regExR_on_slope,3);
out.regExR_on_slope_pd_finiteSample = mean(regExR_on_slope_pd,3);



end